function [ peakAmplitudes peakTimes peakRiseTimes summationFactor response t ] = summarizeAlphaTrains( traceFilePrefix, traceNums, delay, ISI, numPulses, showFigure )
%ANALYZEALPHATRAIN Summarizes peak values of response to train of
%alpha-synapse current stimulation, across multiple responses
%
% Measures the amplitude above baseline of each peak in the response to a
% train of alpha synapses. Assumes positive-going responses. 
% 
% Amplitudes are normalized to the first pulse in the train.
%
% Arguments:
%  traceFile    e.g. 'BS0090/BS0090AAAA'
%  traceNums    e.g. [45 47:48]
%  delay        time to first pulse onset, in sec, e.g. 0.1
%  ISI          time between puls onsets, in sec, e.g. 0.05
%  numPulses    # of alpha-synapse input pulses
% 
% 2009-10-09 Ben Suter
% 2009-10-12 Ben Suter -- display summation factor on figure
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

for f=1:length(traceNums)
    xsg = [ traceFilePrefix sprintf('%04.0f', traceNums(f)) '.xsg' ];
    [ pa peakTimes prt sf r t ] = analyzeAlphaTrain( xsg, delay, ISI, numPulses);
    peakAmplitudes(f,:) = pa;
    peakRiseTimes(f,:) = prt;
    summationFactor(f) = sf;
    response(f,:) = r;
	t(f,:) = t;
end

% % average the repeated responses
% response = mean(response, 1);

% normalize to first peak (in each row)
peakAmplitudes = peakAmplitudes ./ repmat(peakAmplitudes(:,1), 1, size(peakAmplitudes, 2));

% Plot average, normalized peak values with error bars
if showFigure
    figure;
    if length(traceNums) > 1 
        errorbar(peakTimes, mean(peakAmplitudes, 1), std(peakAmplitudes, 1), '.');
    else
        plot(peakTimes, mean(peakAmplitudes, 1), '.');
    end
    yl = ylim(); yl(1) = 0; ylim(yl); % y-axis starts from zero
    xlabel('Time (sec)'); ylabel('Peak response, normalized to first peak (mV)');
    title(sprintf('Response to train of alpha-synapse current injections (N = %g)', length(traceNums)));
    % Display the summation factor next to the last data point
    text(peakTimes(end) + ISI/4, peakAmplitudes(end), {'Summation:', num2str(mean(summationFactor))});
end

end

